args <- commandArgs(TRUE)
setwd(args[1])
#library(genefilter,lib.loc=args[2])
#library(gcrma,lib.loc=args[2])
#library(Biobase,lib.loc=args[2])
#library(preprocessCore,lib.loc=args[2])
#library(MASS,lib.loc=args[2])
#library(affy,lib.loc=args[2])
#library(simpleaffy,lib.loc=args[2])

#library(frma,lib.loc=args[2])
#library(hgu133plus2geneccdscdf)
#library(hgu133plus2frmavecs,lib.loc=args[2])
#library(affyPLM)
.libPaths("/usr/local/biotools/r/R-2.14.0/lib64/R/library")
library(affy)
library(simpleaffy)
library(genefilter)
library(gcrma)
library(Biobase)
library(preprocessCore)
library(MASS)
library(affyPLM)

.libPaths("/home/m081429/R/x86_64-unknown-linux-gnu-library/2.12")
library(frma)
library(hgu133plus2frmavecs)
data(hgu133plus2frmavecs)

frma1 <- function(object, background="rma", normalize="quantile", summarize="robust_weighted_average", target="core", input.vecs=list(normVec=NULL, probeVec=NULL, probeVarBetween=NULL, probeVarWithin=NULL, probesetSD=NULL), output.param=NULL, verbose=FALSE){

  if(!class(object) %in% c("AffyBatch", "ExonFeatureSet")) stop("object must be of class AffyBatch or ExonFeatureSet.")

  if(class(object)=="ExonFeatureSet") target <- match.arg(target, c("core", "full", "extended", "probeset"))
  if(verbose & class(object)=="ExonFeatureSet") message(paste("Exon summarization at ", target, " level.\n", sep=""))

  if(!background %in% c("none", "rma")) stop("background must be either none or rma")
  if(!normalize %in% c("none", "quantile")) stop("normalize must be either none or quantile")
  if(!summarize %in% c("median_polish", "average", "median", "weighted_average", "robust_weighted_average", "random_effect")) stop("summarize must be one of: median_polish, average, median, weighted_average, robust_weighted_average, random_effect")

  if(summarize=="random_effect" & (any(input.vecs$probeVarBetween==0) | any(input.vecs$probeVarWithin==0))) stop("If summarize method is random_effect then probeVarBetween and probeVarWithin must be greater than zero for all probes.")

  if(summarize=="random_effect" & class(object)=="ExonFeatureSet") stop("Summarization method random_effect is not implemented for ExonFeatureSet objects.")
  
  if(class(object)=="AffyBatch"){
    cdfname <- cleancdfname(cdfName(object))
    platform <- gsub("cdf","",cdfname)
  }
  if(class(object)=="ExonFeatureSet"){
    cdfname <- annotation(object)
    platform <- paste(cdfname, target, sep="")
  }

  if(summarize == "median_polish") output <- frmaMedPol(object, background, normalize, target, input.vecs, verbose)
  if(summarize %in% c("average", "median", "weighted_average", "robust_weighted_average")) output <- frmaRobReg1(object, background, normalize, summarize, target, input.vecs, output.param, verbose)
  if(summarize == "random_effect") output <- frmaBatch(object, background, normalize, input.vecs, output.param, verbose)

  if(is.null(output.param)){
    if(is.null(output$stderr)){
      e <- new("ExpressionSet", assayData=assayDataNew(exprs=output$exprs), annotation=platform)
    } else{
      e <- new("ExpressionSet", assayData=assayDataNew(exprs=output$exprs, se.exprs=output$stderr), annotation=platform)
    }
  } else {
    if(is.null(output$stderr)){
      e <- new("frmaExpressionSet", assayData=assayDataNew(exprs=output$exprs), annotation=platform)
    } else{
      e <- new("frmaExpressionSet", assayData=assayDataNew(exprs=output$exprs, se.exprs=output$stderr), annotation=platform)
    }
    if("weights" %in% output.param) w <- output$weights else w <- matrix(nrow=0, ncol=0)
    if("residuals" %in% output.param) r <- output$residuals else r <- matrix(nrow=0, ncol=0)
    if("random_effects" %in% output.param) g <- output$gammas else g <- matrix(nrow=0, ncol=0)
    e@weights <- w
    e@residuals <- r
    e@randomeffects <- g
  }

  return(e)
}

frmaRobReg1 <- function(object, background, normalize, summarize, target, input.vecs, output.param, verbose){

  if(class(object)=="AffyBatch") cdfname <- cleancdfname(cdfName(object))
  if(class(object)=="ExonFeatureSet") cdfname <- annotation(object)
  tmp <- gsub("cdf","",cdfname)
  platform <- gsub("..entrezg", "", tmp)
  if(class(object)=="AffyBatch") vecdataname <- paste(tmp, "frmavecs", sep="")
  if(class(object)=="ExonFeatureSet") vecdataname <- paste(tmp, target, "frmavecs", sep="")
  
  if(background == "rma"){
    if(verbose) message("Background Correcting ...\n")
    if(class(object)=="AffyBatch") object <- bg.correct.rma(object)
    if(class(object)=="ExonFeatureSet") object <- backgroundCorrect(object, verbose=FALSE)
  }

  if(class(object)=="ExonFeatureSet"){
    if(target=="probeset"){
      featureInfo <- getFidProbeset(object)
    }
    if(target=="core"){
      featureInfo <- getFidMetaProbesetCore(object)
    }
    if(target=="full"){
      featureInfo <- getFidMetaProbesetFull(object)
    }
    if(target=="extended"){
      featureInfo <- getFidMetaProbesetExtended(object)
    }

    pmi <- featureInfo[["fid"]]
    pns <- as.character(featureInfo[["fsetid"]])
    pms <- exprs(object)[pmi,, drop=FALSE]
  }
  if(class(object)=="AffyBatch"){
    pms <- pm(object)
    pns <- probeNames(object)
  }
  
  if(is.null(input.vecs$normVec) | is.null(input.vecs$probeVec) | is.null(input.vecs$probeVarWithin) | is.null(input.vecs$probeVarBetween) | (summarize=="robust_weighted_average" & is.null(input.vecs$probesetSD))){
    if(class(object)=="AffyBatch") pkg <- paste(platform, "frmavecs", sep="")
    if(class(object)=="ExonFeatureSet") pkg <- paste(platform, target, "frmavecs", sep="")
    require(pkg, character.only=TRUE, quiet=TRUE) || stop(paste(pkg, "package must be installed first"))
    data(list=eval(vecdataname))

    if(is.null(input.vecs$normVec)) input.vecs$normVec <- get(vecdataname)$normVec
    if(is.null(input.vecs$probeVec)) input.vecs$probeVec <- get(vecdataname)$probeVec
    if(is.null(input.vecs$probeVarWithin)) input.vecs$probeVarWithin <- get(vecdataname)$probeVarWithin
    if(is.null(input.vecs$probeVarBetween)) input.vecs$probeVarBetween <- get(vecdataname)$probeVarBetween
    if(is.null(input.vecs$probesetSD) & summarize=="robust_weighted_average") input.vecs$probesetSD <- get(vecdataname)$probesetSD
  }

  if(normalize == "quantile"){
    if(verbose) message("Normalizing ...\n")
    pms <- normalize.quantiles.use.target(pms, input.vecs$normVec)
	kk<-log2(pms)
    write.table(kk,"afternormalization",sep="\t",row.names=FALSE,col.names=FALSE);
  
	}

  if(verbose) message("Summarizing ...\n")
  pms <- log2(pms)
  
  if(summarize == "average"){  
    exprs <- subColSummarizeAvg((pms-input.vecs$probeVec), pns)
    weights <- NULL
    stderr <- NULL
  }

  if(summarize == "median"){  
    exprs <- subColSummarizeMedian((pms-input.vecs$probeVec), pns)
    weights <- NULL
    stderr <- NULL
  }
  
  if(summarize == "weighted_average"){  
    w <- 1/(input.vecs$probeVarWithin + input.vecs$probeVarBetween)
    if(any(input.vecs$probeVarWithin==0) | any(input.vecs$probeVarBetween==0)) message("Either probeVarWithin or probeVarBetween is 0 for some probes -- setting corresponding weights to 1")
    w[w==Inf] <- 1
    exprs <- subColSummarizeAvg((pms-input.vecs$probeVec)*w, pns)
    W <- as.vector(rowsum(w, pns, reorder=FALSE))
    exprs <- (exprs/W)*as.vector(rowsum(rep(1,length(pns)),pns,reorder=FALSE))
    weights <- NULL
    stderr <- NULL
  }

  if(summarize == "robust_weighted_average"){
    w <- 1/(input.vecs$probeVarWithin + input.vecs$probeVarBetween)
    if(any(input.vecs$probeVarWithin==0) | any(input.vecs$probeVarBetween==0)) message("Either probeVarWithin or probeVarBetween is 0 for some probes -- setting corresponding weights to 1")
    w[w==Inf] <- 1
    N <- 1:dim(pms)[1]
    S <- split(N, pns)
    fit <- lapply(1:length(S), function(i) {
	s <- S[[i]]
	rwaFit2(pms[s,, drop=FALSE], w[s], input.vecs$probeVec[s], input.vecs$probesetSD[i])
    })
    names(fit) <- unique(pns)
    exprs <- matrix(unlist(lapply(fit, function(x) x$Estimates)), ncol=ncol(pms), byrow=TRUE)
    rownames(exprs) <- names(fit)
    colnames(exprs) <- colnames(pms)
    if("weights" %in% output.param){
      weights <- matrix(unlist(lapply(fit, function(x) x$Weights)), ncol=ncol(pms), byrow=TRUE)
      rownames(weights) <- pns
      colnames(weights) <- sampleNames(object)
    } else weights <- NULL
    stderr <- matrix(unlist(lapply(fit, function(x) x$StdErrors)), ncol=ncol(pms), byrow=TRUE)
    rownames(stderr) <- names(fit)
    colnames(stderr) <- sampleNames(object)
  }

  if("residuals" %in% output.param){
    residuals <- apply(exprs,2,function(x) rep(x, table(pns)))
    residuals <- (pms-input.vecs$probeVec) - residuals
    rownames(residuals) <- pns
    colnames(residuals) <- sampleNames(object)
  } else residuals <- NULL
  
  colnames(exprs) <- sampleNames(object)
  
  return(list(exprs=exprs, stderr=stderr, weights=weights, residuals=residuals, gammas=NULL))
}

rwaFit2 <- function(x1, x2, x3, x4){
  ncols <- ncol(x1)
  w.tmp <- x2/max(x2)
  w.tmp <- matrix(rep(w.tmp, ncols), ncol=ncols)
  pe.tmp <- x3
  pe.tmp[1] <- pe.tmp[1]-sum(pe.tmp)
  rcModelWPLM(y=x1, w=w.tmp, row.effects=pe.tmp, input.scale=x4)
}


data<-ReadAffy()
before<-log2(pm(data))
d<-frma1(data, normalize = "quantile",summarize = "average")
after<-read.table("afternormalization",sep="\t")
diff<-after-before
meandiff<-mean(diff)
all<-abs(diff[,1]-meandiff[1])
length(all)
i<-2
len<-dim(diff)
while(i<len[2]+1)
{
        madvec<- abs(diff[,i]-meandiff[i])
        all<-cbind(all,madvec)
        i=i+1
}
finalmean<-apply(all,2,mean)
col<-colnames(before)
cbind(col,finalmean)
fin<-as.matrix(cbind(col,finalmean))
write.table(fin,"result.txt",sep="\t",quote=FALSE,row.names =FALSE,col.names =FALSE)

